library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
df <- read.csv("../data/study1.csv") |>
mutate(
Date = as.Date(Date),
Participant = fct_reorder(Participant, Date),
Screen_Refresh = as.character(Screen_Refresh),
Illusion_Side = as.factor(Illusion_Side),
Block = as.factor(Block),
Education = fct_relevel(Education, "Master", "Bachelor", "High School", "Other")
)
outliers <- c(
# Half of the trials have very short RT
# Prolific Status: REJECTED
"S33",
# Block n2 with very short RTs
# Prolific Status: RETURNED
"S20",
# Error rate of 46.2% and short RTs
# Prolific Status: RETURNED
"S51",
# Error rate of 46.2%
# Prolific Status: RETURNED
"S49",
# Error rate of 47.9%
# Prolific Status: REJECTED
"S47",
# Error rate of 42.1% and very large RT SD
# Prolific Status: REJECTED
"S12"
)
We removed 6 participants upon inspection of the average error rage (when close to 50%, suggesting random answers) and/or when the reaction time distribution was implausibly fast.
For each block, we computed the error rate and, if more than 50%, we discarded the whole block (as it likely indicates that instructions got mixed up, for instance participants were selecting the smaller instead of the bigger circle).
dfsub <- df |>
group_by(Participant) |>
summarize(
# n = n(),
Error = sum(Error) / n(),
RT_Mean = mean(RT),
RT_SD = sd(RT),
) |>
ungroup() |>
arrange(desc(Error))
knitr::kable(dfsub) |>
kableExtra::row_spec(which(dfsub$Participant %in% outliers), background = "#EF9A9A")
| Participant | Error | RT_Mean | RT_SD |
|---|---|---|---|
| S47 | 0.479 | 262 | 296 |
| S51 | 0.465 | 263 | 372 |
| S49 | 0.462 | 630 | 611 |
| S12 | 0.421 | 507 | 1679 |
| S20 | 0.402 | 492 | 725 |
| S41 | 0.300 | 723 | 579 |
| S50 | 0.287 | 659 | 333 |
| S16 | 0.269 | 578 | 172 |
| S05 | 0.262 | 716 | 271 |
| S04 | 0.262 | 588 | 206 |
| S42 | 0.260 | 718 | 1294 |
| S43 | 0.255 | 837 | 711 |
| S33 | 0.251 | 599 | 1326 |
| S46 | 0.246 | 699 | 414 |
| S06 | 0.246 | 560 | 197 |
| S39 | 0.243 | 693 | 232 |
| S15 | 0.238 | 1160 | 1660 |
| S40 | 0.227 | 748 | 243 |
| S08 | 0.225 | 780 | 427 |
| S27 | 0.219 | 785 | 836 |
| S30 | 0.216 | 506 | 150 |
| S26 | 0.215 | 738 | 375 |
| S10 | 0.213 | 1076 | 557 |
| S23 | 0.210 | 870 | 489 |
| S52 | 0.207 | 804 | 378 |
| S34 | 0.206 | 652 | 367 |
| S28 | 0.205 | 511 | 147 |
| S32 | 0.203 | 627 | 223 |
| S35 | 0.202 | 1133 | 982 |
| S44 | 0.202 | 869 | 424 |
| S24 | 0.201 | 720 | 298 |
| S01 | 0.198 | 1110 | 738 |
| S48 | 0.196 | 867 | 652 |
| S18 | 0.189 | 983 | 830 |
| S45 | 0.187 | 711 | 195 |
| S25 | 0.186 | 803 | 397 |
| S36 | 0.185 | 846 | 765 |
| S09 | 0.182 | 1045 | 1082 |
| S19 | 0.181 | 1001 | 562 |
| S22 | 0.176 | 1307 | 1091 |
| S14 | 0.175 | 1146 | 900 |
| S02 | 0.174 | 703 | 243 |
| S29 | 0.173 | 1046 | 918 |
| S38 | 0.173 | 616 | 381 |
| S17 | 0.171 | 927 | 550 |
| S37 | 0.170 | 769 | 312 |
| S13 | 0.157 | 848 | 364 |
| S31 | 0.156 | 712 | 383 |
| S07 | 0.142 | 1109 | 863 |
| S11 | 0.139 | 895 | 816 |
| S03 | 0.124 | 741 | 331 |
| S21 | 0.092 | 884 | 509 |
temp <- df |>
group_by(Participant, Illusion_Type, Block) |>
summarize(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
arrange(desc(ErrorRate_per_block))
temp2 <- temp |>
filter(ErrorRate_per_block >= 0.5) |>
group_by(Illusion_Type, Block) |>
summarize(n = n()) |>
arrange(desc(n), Illusion_Type) |>
ungroup() |>
mutate(n_trials = cumsum(n * 56),
p_trials = n_trials / nrow(df))
# knitr::kable(temp2)
p1 <- temp |>
estimate_density(at = c("Illusion_Type", "Block")) |>
ggplot(aes(x = x, y = y)) +
geom_line(aes(color = Illusion_Type, linetype = Block)) +
geom_vline(xintercept = 0.5, linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
labs(y = "Distribution", x = "Error Rate") +
theme_modern()
p2 <- temp2 |>
mutate(Block = fct_rev(Block)) |>
ggplot(aes(x = Illusion_Type, y = p_trials)) +
geom_bar(stat="identity", aes(fill = Block)) +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(labels = scales::percent, expand = c(0, 0)) +
labs(y = "Percentage of Trials Removed", x = "Illusion Type") +
theme_modern() +
theme(axis.text.x = element_text(angle=45, hjust = 1))
p1 | p2
# Drop
df <- df |>
group_by(Participant, Illusion_Type, Block) |>
mutate(ErrorRate_per_block = sum(Error) / n()) |>
ungroup() |>
filter(ErrorRate_per_block < 0.5) |>
select(-ErrorRate_per_block)
rm(temp, temp2)
# RT distribution
p <- estimate_density(df, select = "RT", at = c("Participant", "Block")) |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
mutate(color = ifelse(Participant %in% outliers, "red", "blue")) |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = color, group = interaction(Participant, Block), linetype = Block)) +
geom_vline(xintercept = 2500, linetype = "dashed", color = "red") +
scale_color_manual(values=c("red"="red", "blue"="blue"), guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(xlim = c(0, 3000)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
facet_wrap(~Participant) +
labs(y = "", x = "Reaction Time (ms)")
p
# ggsave("figures/outliers.png", p, width=20, height=15)
# Filter out
df <- filter(df, !Participant %in% outliers)
p1 <- estimate_density(df, select = "RT", at = "Participant") |>
group_by(Participant) |>
normalize(select = "y") |>
ungroup() |>
ggplot(aes(x = x, y = y)) +
geom_area(data = normalize(estimate_density(df, select = "RT"), select = "y"), alpha = 0.2) +
geom_line(aes(color = Participant, group = Participant)) +
geom_vline(xintercept = c(150, 3000), linetype = "dashed", color = "red") +
scale_color_material_d("rainbow", guide = "none") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
coord_cartesian(xlim = c(0, 3500)) +
theme_modern() +
theme(axis.text.y = element_blank()) +
# facet_wrap(~Participant) +
labs(y = "", x = "Reaction Time (ms)")
df$Outlier <- df$RT < 150 | df$RT > 3000
p2 <- df |>
group_by(Participant) |>
summarize(Outlier = sum(Outlier) / n()) |>
mutate(Participant = fct_reorder(Participant, Outlier)) |>
ggplot(aes(x = Participant, y = Outlier)) +
geom_bar(stat = "identity", aes(fill = Participant)) +
scale_fill_material_d("rainbow", guide = "none") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0), labels = scales::percent) +
see::theme_modern() +
theme(axis.text.x = element_blank())
p1 | p2
We removed 692 (1.37%) outlier trials (150 ms < RT < 3000 ms).
df <- filter(df, Outlier == FALSE)
dfsub <- df |>
group_by(Participant) |>
select(Participant, Age, Sex, Education, Nationality, Ethnicity, Duration, Break_Duration, Screen_Resolution, Screen_Refresh, Device_OS) |>
slice(1) |>
ungroup()
The final sample included 46 participants (Mean age = 26.7, SD = 7.7, range: [19, 60]; Sex: 39.1% females, 56.5% males, 4.3% other; Education: Master, 17.39%; Bachelor, 41.30%; High School, 39.13%; Other, 2.17%).
plot_distribution <- function(dfsub, what = "Age", title = what, subtitle = "", fill = "orange") {
dfsub |>
ggplot(aes_string(x = what)) +
geom_density(fill = fill) +
geom_vline(xintercept = mean(dfsub[[what]]), color = "red", linetype = "dashed") +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
ggtitle(title, subtitle = subtitle) +
theme_modern() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(face = "italic", hjust = 0.5),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank()
)
}
plot_waffle <- function(dfsub, what = "Nationality") {
ggwaffle::waffle_iron(dfsub, what) |>
# mutate(label = emojifont::fontawesome('fa-twitter')) |>
ggplot(aes(x, y, fill = group)) +
ggwaffle::geom_waffle() +
# geom_point() +
# geom_text(aes(label=label), family='fontawesome-webfont', size=4) +
coord_equal() +
ggtitle(what) +
labs(fill = "") +
theme_void() +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
p1 <- plot_distribution(dfsub, "Age", fill = "#FF9800")
p2 <- plot_distribution(dfsub, "Duration", title = "Total Duration", subtitle = "in minutes", fill = "#F44336")
p3 <- plot_distribution(dfsub, "Break_Duration", title = "Break Duration", subtitle = "in minutes", fill = "#3F51B5")
p4 <- plot_waffle(dfsub, "Sex") +
scale_fill_manual(values = c("Male" = "#2196F3", "Female" = "#E91E63", "Other" = "#FF9800"))
p5 <- plot_waffle(dfsub, "Education") +
scale_fill_viridis_d()
p6 <- plot_waffle(dfsub, "Nationality") +
scale_fill_metro_d()
p7 <- plot_waffle(dfsub, "Ethnicity") +
scale_fill_manual(values = c("Latino" = "#FF5722", "Asian" = "#FF9800", "Caucasian" = "#2196F3", "African" = "#4CAF50", "Jewish" = "#9C27B0"))
p8 <- plot_waffle(dfsub, "Screen_Resolution") +
scale_fill_pizza_d()
p9 <- plot_waffle(dfsub, "Device_OS") +
scale_fill_bluebrown_d()
# p10 <- plot_waffle(dfsub, "Screen_Refresh") +
# scale_fill_viridis_d()
(p1 / p2 / p3) | (p4 / p5 / p6) | (p7 / p8 / p9)
plot_descriptive <- function(data, side="leftright") {
if(side == "leftright") {
x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
x <- tools::toTitleCase(gsub("arrow", "", x))
if(x == "Left") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Left", "Right")
} else if(x == "Right") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Right", "Left")
}
} else {
x <- data[data$Error == 0 & data$Illusion_Side == 1, ]$Answer[1]
x <- tools::toTitleCase(gsub("arrow", "", x))
if(x == "Up") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Up", "Down")
} else if(x == "Down") {
data$Answer <- ifelse(data$Illusion_Side == 1, "Down", "Up")
}
data$Answer <- fct_rev(data$Answer)
}
dodge1 <- 0.1 * diff(range(data$Illusion_Difference))
dodge2 <- -0.1 * diff(range(data$Illusion_Strength))
colors <- colorRampPalette(c("#4CAF50", "#009688", "#00BCD4", "#2196F3", "#3F51B5", "#673AB7", "#9C27B0"))(length(unique(data$Illusion_Strength)))
p1 <- data |>
group_by(Illusion_Difference, Illusion_Strength, Answer) |>
summarize(Error = mean(Error)) |>
mutate(Illusion_Strength = as.factor(round(Illusion_Strength, 2))) |>
ggplot(aes(x = Illusion_Difference, y = Error)) +
geom_bar(aes(fill=Illusion_Strength), position = position_dodge(width=dodge1), stat="identity") +
geom_line(aes(color = Illusion_Strength), position = position_dodge(width=dodge1)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
theme_modern() +
labs(
color = "Illusion Strength",
fill = "Illusion Strength",
y = "Probability of Error",
x = "Task Difficulty"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data$Illusion_Difference)))
p2 <- data |>
group_by(Illusion_Difference, Illusion_Strength, Answer) |>
summarize(Error = mean(Error)) |>
mutate(Illusion_Difference = as.factor(round(Illusion_Difference, 2))) |>
ggplot(aes(x = Illusion_Strength, y = Error)) +
geom_vline(xintercept=0, linetype="dotted", alpha=0.6) +
geom_bar(aes(fill=Illusion_Difference), position = position_dodge(width=dodge2), stat="identity") +
geom_line(aes(color = Illusion_Difference), position = position_dodge(width=dodge2)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
theme_modern() +
labs(
color = "Task Difficulty",
fill = "Task Difficulty",
y = "Probability of Error",
x = "Illusion Strength"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
if(side == "leftright") {
p <- ((p1 + facet_wrap(~Answer, ncol=2, labeller = "label_both")) /
(p2 + facet_wrap(~Answer, ncol=2, labeller = "label_both"))) +
plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"),
theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
} else {
p <- ((p1 + facet_wrap(~Answer, nrow=2, labeller = "label_both")) |
(p2 + facet_wrap(~Answer, nrow=2, labeller = "label_both"))) +
plot_annotation(title = paste(data$Illusion_Type[1], "Illusion"),
theme = theme(plot.title = element_text(face = "bold", hjust = 0.5)))
}
p
}
data <- filter(df, Illusion_Type == "Delboeuf")
plot_descriptive(data)
best_models <- function(data) {
models <- list()
for(i in 1:1) {
for(j in 1:1) {
for(k1 in c("", "_log", "_sqrt", "_cbrt")) {
for(k2 in c("")) {
for(side in c("", "-side")) {
name <- paste0("dif", k1, i, "-", "str", k2, j, side)
# print(name)
f <- paste0("poly(Illusion_Difference",
k1,
", ",
i,
") * poly(Illusion_Strength",
k2,
", ",
j,
") + (1|Participant)")
if(side == "-side") f <- paste0("Illusion_Side * ", f)
m <- glmmTMB::glmmTMB(as.formula(paste0("Error ~ ", f)),
data = data, family = "binomial")
if(performance::check_convergence(m)) {
models[[name]] <- m
}
}
}
}
}
}
to_keep <- compare_performance(models, metrics = c("BIC")) |>
arrange(BIC) |>
slice(1:10) |>
pull(Name)
test <- test_performance(models[to_keep], reference=1)
perf <- compare_performance(models[to_keep], metrics = c("BIC", "R2"))
merge(perf, test) |>
arrange(BIC) |>
select(Name, BIC, R2_marginal, BF) |>
mutate(BF = insight::format_bf(BF, name=""))
}
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 3305 0.392
## 2 dif_sqrt1-str1 3307 0.400 = 0.353
## 3 dif_log1-str1 3317 0.412 = 0.003
## 4 dif_cbrt1-str1-side 3330 0.394 < 0.001
## 5 dif_sqrt1-str1-side 3332 0.401 < 0.001
## 6 dif1-str1 3333 0.424 < 0.001
## 7 dif_log1-str1-side 3341 0.414 < 0.001
## 8 dif1-str1-side 3356 0.426 < 0.001
cbrt <- function(x) sign(x) * abs(x)**(1/3)
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 12.5 seconds.
## Chain 3 finished in 14.5 seconds.
## Chain 1 finished in 15.4 seconds.
## Chain 2 finished in 16.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 14.6 seconds.
## Total execution time: 16.2 seconds.
# parameters::parameters(model)
plot_model <- function(data, model) {
data <- mutate(data, .dots_side = ifelse(Error == 1, "bottom", "top"))
# Get variables
vars <- insight::find_predictors(model)$conditional
vardiff <- vars[1]
varstrength <- vars[2]
# Get predicted
pred <- estimate_relation(model,
at = vars,
length = c(NA, 25))
pred[[vardiff]] <- as.factor(pred[[vardiff]])
# Set colors for lines
colors <- colorRampPalette(c("#F44336", "#FFC107", "#4CAF50"))(length(unique(data[[vardiff]])))
diffvals <- as.numeric(as.character(unique(sort(pred[[vardiff]]))))
names(colors) <- diffvals
# Assign color from the same palette to every observation of data (for geom_dots)
closest <- diffvals[max.col(-abs(outer(data[[vars[1]]], diffvals, "-")))]
data$color <- colors[as.character(closest)]
data$color <- fct_reorder(data$color, closest)
# Manual jittering
xrange <- 0.05*diff(range(data[[varstrength]]))
data$x <- data[[varstrength]]
data$x[data$x > 0] <- data$x[data$x > 0] - runif(sum(data$x > 0), 0, xrange)
data$x[data$x < 0] <- data$x[data$x < 0] + runif(sum(data$x < 0), 0, xrange)
data$x[round(data$x, 2) == 0] <- data$x[round(data$x, 2) == 0] + runif(sum(round(data$x, 2) == 0), -xrange/2, xrange/2)
pred |>
ggplot(aes_string(x = varstrength, y = "Predicted")) +
geom_dots(
data = data,
aes(x=x, y = Error, group = Error, side = .dots_side, order=color),
fill = data$color,
color = NA,
alpha=0.5) +
geom_slab(data=data, aes(y = Error)) +
geom_ribbon(aes_string(ymin = "CI_low", ymax = "CI_high", fill = vardiff, group = vardiff), alpha = 0.2) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_hline(yintercept = c(0.05, 0.5, 0.95), linetype = "dotted", alpha = 0.5) +
geom_line(aes_string(color = vardiff, group = vardiff)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0), labels = scales::percent) +
scale_x_continuous(expand = c(0, 0)) +
scale_color_manual(values = colors) +
scale_fill_manual(values = colors) +
coord_cartesian(xlim=c(min(data[[varstrength]]), max(data[[varstrength]]))) +
theme_modern() +
labs(
title = paste0(data$Illusion_Type[1], " Illusion"),
color = "Difficulty", fill = "Difficulty",
y = "Probability of Error",
x = "Illusion Strength"
) +
theme(plot.title = element_text(face = "bold", hjust = 0.5))
}
plot_model(data, model)
make_gam <- function(data) {
formula <- brms::bf(
Error ~ t2(Illusion_Difference, Illusion_Strength, bs = "cr", k=4) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
list(p = plot_model(data, model), model = model)
}
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 28.2 seconds.
## Chain 3 finished in 28.5 seconds.
## Chain 1 finished in 28.9 seconds.
## Chain 2 finished in 30.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 28.9 seconds.
## Total execution time: 30.2 seconds.
gam$p
# Parameter Standardization
std_params <- function(model, min=0, max=2) {
estimate_relation(
model,
at = list(Illusion_Strength = c(0),
Illusion_Difference = seq(min, max, length.out=500)),
) |>
select(Illusion_Strength, Illusion_Difference, Error = Predicted) |>
slice(c(which.min(abs(Error - 0.005)),
which.min(abs(Error - 0.025)),
which.min(abs(Error - 0.25)))) |>
mutate(Error = insight::format_value(Error, as_percent=TRUE))
}
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=2)
std_params(gam$model, min=0, max=2)
data <- filter(df, Illusion_Type == "Ebbinghaus")
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2212 0.615
## 2 dif_sqrt1-str1 2215 0.615 = 0.199
## 3 dif_cbrt1-str1-side 2227 0.623 < 0.001
## 4 dif_sqrt1-str1-side 2230 0.623 < 0.001
## 5 dif_log1-str1 2232 0.617 < 0.001
## 6 dif1-str1 2258 0.619 < 0.001
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 18.4 seconds.
## Chain 3 finished in 18.6 seconds.
## Chain 2 finished in 18.8 seconds.
## Chain 1 finished in 19.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 18.7 seconds.
## Total execution time: 19.2 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 2 finished in 29.6 seconds.
## Chain 4 finished in 30.3 seconds.
## Chain 1 finished in 32.0 seconds.
## Chain 3 finished in 33.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 31.3 seconds.
## Total execution time: 33.2 seconds.
gam$p
# Parameter Standardization
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=2)
std_params(gam$model, min=0, max=2)
data <- filter(df, Illusion_Type == "Rod-Frame")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 15)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2070 0.464
## 2 dif_log1-str1 2072 0.458 = 0.371
## 3 dif_cbrt1-str1 2075 0.452 = 0.063
## 4 dif_sqrt1-str1-side 2095 0.478 < 0.001
## 5 dif1-str1 2095 0.505 < 0.001
## 6 dif_log1-str1-side 2097 0.471 < 0.001
## 7 dif_cbrt1-str1-side 2101 0.465 < 0.001
## 8 dif1-str1-side 2119 0.525 < 0.001
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 6.4 seconds.
## Chain 3 finished in 6.6 seconds.
## Chain 4 finished in 6.6 seconds.
## Chain 2 finished in 6.7 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 6.5 seconds.
## Total execution time: 6.9 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 13.8 seconds.
## Chain 2 finished in 16.5 seconds.
## Chain 4 finished in 17.1 seconds.
## Chain 3 finished in 17.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 16.2 seconds.
## Total execution time: 17.5 seconds.
gam$p
# Parameter Standardization
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0.01, max=12)
std_params(gam$model, min=0.01, max=12)
data <- filter(df, Illusion_Type == "Vertical-Horizontal")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 90)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2011 0.661
## 2 dif_cbrt1-str1 2014 0.657 = 0.234
## 3 dif_log1-str1 2017 0.674 = 0.050
## 4 dif1-str1 2022 0.679 = 0.004
## 5 dif_sqrt1-str1-side 2031 0.665 < 0.001
## 6 dif_cbrt1-str1-side 2034 0.662 < 0.001
## 7 dif_log1-str1-side 2038 0.677 < 0.001
## 8 dif1-str1-side 2043 0.681 < 0.001
formula <- brms::bf(
Error ~ sqrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 22.4 seconds.
## Chain 1 finished in 22.6 seconds.
## Chain 3 finished in 22.8 seconds.
## Chain 2 finished in 23.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 22.9 seconds.
## Total execution time: 24.0 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 3 finished in 22.6 seconds.
## Chain 1 finished in 22.7 seconds.
## Chain 4 finished in 23.1 seconds.
## Chain 2 finished in 23.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 23.0 seconds.
## Total execution time: 23.8 seconds.
gam$p
# Parameter Standardization
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.40)
std_params(gam$model, min=0, max=0.40)
data <- filter(df, Illusion_Type == "Zöllner")
plot_descriptive(data)
data <- filter(data, abs(Illusion_Strength) < 45)
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 1041 0.405
## 2 dif_log1-str1 1043 0.414 = 0.464
## 3 dif_sqrt1-str1 1044 0.423 = 0.260
## 4 dif1-str1 1065 0.491 < 0.001
## 5 dif_cbrt1-str1-side 1066 0.409 < 0.001
## 6 dif_log1-str1-side 1068 0.418 < 0.001
## 7 dif_sqrt1-str1-side 1069 0.427 < 0.001
## 8 dif1-str1-side 1090 0.494 < 0.001
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * Illusion_Strength +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 9.8 seconds.
## Chain 4 finished in 10.1 seconds.
## Chain 2 finished in 10.7 seconds.
## Chain 3 finished in 12.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 10.6 seconds.
## Total execution time: 12.0 seconds.
# parameters::parameters(model)
plot_model(data, model)
# Parameter Standardization
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=5)
data <- filter(df, Illusion_Type == "White")
plot_descriptive(data)
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_log1-str1 2176 0.783
## 2 dif_cbrt1-str1 2178 0.784 = 0.404
## 3 dif_sqrt1-str1 2181 0.784 = 0.101
## 4 dif_log1-str1-side 2188 0.788 = 0.003
## 5 dif_cbrt1-str1-side 2190 0.789 = 0.001
## 6 dif_sqrt1-str1-side 2192 0.790 < 0.001
## 7 dif1-str1 2195 0.787 < 0.001
## 8 dif1-str1-side 2206 0.793 < 0.001
formula <- brms::bf(
Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 79.3 seconds.
## Chain 2 finished in 89.6 seconds.
## Chain 3 finished in 90.0 seconds.
## Chain 4 finished in 93.3 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 88.0 seconds.
## Total execution time: 93.5 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 3 finished in 39.7 seconds.
## Chain 1 finished in 43.7 seconds.
## Chain 4 finished in 51.6 seconds.
## Chain 2 finished in 52.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 46.8 seconds.
## Total execution time: 52.2 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
# unique(data$Illusion_Strength)
std_params(model, min=0, max=20)
std_params(gam$model, min=0, max=20)
data <- filter(df, Illusion_Type == "Müller-Lyer")
plot_descriptive(data, side = "updown")
best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2449 0.749
## 2 dif_log1-str1-side 2534 0.748 < 0.001
## 3 dif_cbrt1-str1-side 2535 0.757 < 0.001
formula <- brms::bf(
Error ~ log(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 2 finished in 16.7 seconds.
## Chain 4 finished in 18.4 seconds.
## Chain 3 finished in 19.4 seconds.
## Chain 1 finished in 19.6 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 18.5 seconds.
## Total execution time: 19.7 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 32.9 seconds.
## Chain 2 finished in 33.8 seconds.
## Chain 3 finished in 33.8 seconds.
## Chain 1 finished in 34.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 33.6 seconds.
## Total execution time: 34.1 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.6)
std_params(gam$model, min=0, max=0.6)
data <- filter(df, Illusion_Type == "Ponzo")
plot_descriptive(data, side = "updown")
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 2551 0.661
## 2 dif_sqrt1-str1 2555 0.665 = 0.187
## 3 dif_log1-str1 2576 0.678 < 0.001
## 4 dif1-str1 2590 0.684 < 0.001
## 5 dif_log1-str1-side 2592 0.687 < 0.001
## 6 dif1-str1-side 2606 0.694 < 0.001
formula <- brms::bf(
Error ~ poly(Illusion_Difference, 2) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 12.3 seconds.
## Chain 2 finished in 13.2 seconds.
## Chain 1 finished in 13.3 seconds.
## Chain 3 finished in 14.7 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 13.4 seconds.
## Total execution time: 14.9 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 29.3 seconds.
## Chain 2 finished in 30.4 seconds.
## Chain 3 finished in 30.3 seconds.
## Chain 4 finished in 35.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 31.3 seconds.
## Total execution time: 35.3 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=20)
std_params(gam$model, min=0, max=20)
data <- filter(df, Illusion_Type == "Poggendorff")
plot_descriptive(data, side = "updown")
data <- filter(data, abs(Illusion_Strength) < 45)
plot_descriptive(data, side = "updown")
best_models(data)
## Name BIC R2_marginal BF
## 1 dif_cbrt1-str1 1900 0.527
## 2 dif_sqrt1-str1 1905 0.526 = 0.086
## 3 dif_log1-str1 1923 0.522 < 0.001
## 4 dif_cbrt1-str1-side 1929 0.529 < 0.001
## 5 dif1-str1 1930 0.521 < 0.001
## 6 dif_sqrt1-str1-side 1934 0.529 < 0.001
## 7 dif_log1-str1-side 1952 0.527 < 0.001
## 8 dif1-str1-side 1958 0.526 < 0.001
formula <- brms::bf(
Error ~ cbrt(Illusion_Difference) * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 2 finished in 11.2 seconds.
## Chain 1 finished in 11.6 seconds.
## Chain 3 finished in 11.8 seconds.
## Chain 4 finished in 11.9 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 11.6 seconds.
## Total execution time: 12.1 seconds.
# parameters::parameters(model)
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 1 finished in 18.1 seconds.
## Chain 4 finished in 19.0 seconds.
## Chain 3 finished in 19.4 seconds.
## Chain 2 finished in 21.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 19.4 seconds.
## Total execution time: 21.1 seconds.
gam$p
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=0.5)
std_params(gam$model, min=0, max=0.5)
data <- filter(df, Illusion_Type == "Contrast")
plot_descriptive(data, side = "updown")
best_models(data)
## Random effect variances not available. Returned R2 does not account for random effects.
## Random effect variances not available. Returned R2 does not account for random effects.
## Name BIC R2_marginal BF
## 1 dif_sqrt1-str1 2353 0.673
## 2 dif_cbrt1-str1 2354 0.673 = 0.713
## 3 dif_log1-str1 2357 0.672 = 0.128
## 4 dif1-str1 2359 0.675 = 0.047
## 5 dif_cbrt1-str1-side 2580 0.685 < 0.001
## 6 dif1-str1-side 2585 0.685 < 0.001
formula <- brms::bf(
Error ~ Illusion_Difference * poly(Illusion_Strength, 2) +
(1 | Participant),
family = "bernoulli"
)
model <- brms::brm(formula,
data = data,
refresh = 0
)
## Running MCMC with 4 parallel chains...
##
## Chain 4 finished in 22.7 seconds.
## Chain 3 finished in 23.0 seconds.
## Chain 1 finished in 23.3 seconds.
## Chain 2 finished in 24.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 23.4 seconds.
## Total execution time: 24.5 seconds.
plot_model(data, model)
gam <- make_gam(data)
## Running MCMC with 4 parallel chains...
##
## Chain 2 finished in 30.7 seconds.
## Chain 1 finished in 31.8 seconds.
## Chain 3 finished in 33.8 seconds.
## Chain 4 finished in 34.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 32.7 seconds.
## Total execution time: 34.6 seconds.
gam$p
# Parameter Standardization
# range(data$Illusion_Difference)
# range(data$Illusion_Strength)
std_params(model, min=0, max=25)
std_params(gam$model, min=0, max=25)